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Abstract 

We are interested in the asymptotic behavior of orthogonal poly¬ 
nomials of the generalized Jacobi type as their degree n goes to oo. 
These are defined on the interval [—1,1] with weight function 

w{x) = {1 — x)°'{l + x)^h{x), a,(3>—l 

and h{x) a real, analytic and strictly positive function on [—1,1]. This 
information is available in the work of Kuijlaars, McLaughlin, Van 
Assche and Vanlessen [22], where the authors use the Riemann-Hilbert 
formulation and the Deift-Zhou non-linear steepest descent method. 

We show that computing higher-order terms can be simplified, lead¬ 
ing to their efficient construction. The resulting asymptotic expansions 
in every region of the complex plane are implemented both symboli¬ 
cally and numerically, and the code is made publicly available. 
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The main advantage of these expansions is that they lead to in¬ 
creasing accuracy for increasing degree of the polynomials, at a com¬ 
putational cost that is actually independent of the degree. In con¬ 
trast, the typical use of the recurrence relation for orthogonal poly¬ 
nomials in computations leads to a cost that is at least linear in the 
degree. Furthermore, the expansions may be used to compute Gaus¬ 
sian quadrature rules in 0{n) operations, rather than 0{ri^) based on 
the recurrence relation. 


1 Introduction 


In this paper we are interested in the symbolic implementation and numeri¬ 
cal computation of asymptotic expansions for monic polynomials 7Tn(x) that 
are orthogonal with respect to a Jacobi-type weight function on the interval 
[- 1 , 1 ]: 




with 


TTn{x)TTk{x)w{x)dx = 0, /c = 0,1,..., n — 1, 
w{x) = (1 — x)'^{l + x)^h{x), a,l3 > —1, 


( 1 . 1 ) 


and where h{x) is a real analytic and strictly positive function on [—1,1]. 

The work of Kuijlaars, McLaughlin, Van Assche and Vanlessen [22] pro¬ 
vides complete asymptotic information of the large n behavior of 7r„(x), 
together with associated quantities such as recurrence coefficients On and 
l^n of the three term recurrence relation 


vrn+l(x) = (x - an)7Tn{x) - /3n7rn-l{x), (1.2) 

as well as leading term coefficients and Hankel determinants. These results 
are obtained using the Riemann-Hilbert formulation for 7rn(x), see the semi¬ 
nal paper of Fokas, Its and Kitaev, [8], and the steepest descent method due 
to Deift and Zhou, [5, 6]. This procedure gives three types of asymptotic 
expansions: outer asymptotics, valid for x G C \ [—1,1], inner asymptotics, 
for X G (—1,1) but away from the endpoints, and boundary asymptotics, for 
|x =F 1 ] < d, for some fixed 6 > 0. 

The purpose of this paper is to provide an automatic and efficient im¬ 
plementation (symbolic and numerical) of the asymptotic expansions for 
these Jacobi-type polynomials 7r„(x) presented in [22]. In this reference, 
the leading order terms are given explicitly, and we detail the derivation of 
higher-order terms. It only requires elementary numerical techniques: in 
particular, it does not need any evaluation of special functions. We will 
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also deal with some specific and non~trivial issues that arise when imple¬ 
menting the formulation of [22] in the complex plane in a symbolic and 
numerical setting. All formulas are implemented in Maple, Julia and 
Matlab and are available on the software webpage of our research group: 
http://nines.cs.kuleuven.be/software/JACOBI. 

From a computational point of view, the asymptotic expansions given 
in [22] present two important advantages; they become increasingly accu¬ 
rate as the degree n becomes large, and the computing time is essentially 
independent of n. In this sense, they compare favourably to other meth¬ 
ods to compute 7r„(a:), such as the use of the recurrence relation (1.2). For 
this reason, the idea of using asymptotic expansions for the computation of 
orthogonal polynomials has already been present in the literature for some 
time, either using approximations coming from integral representations or 
differential equations (see the work of Hale and Townsend [15, 16] and ref¬ 
erences therein), or more recently from Riemann-Hilbert problems (see for 
instance the work of Olver and Trogdon [28, 27, 29, 25, 34]. Additionally, 
these expansions can also be used to construct Gaussian quadrature rules 
with a high number of points: we refer the reader to [13, 3, 15, 30, 2] and 
§ 2 . 6 . 

Observe that when h{x) = 1 we have the standard Jacobi polynomials, 
whose strong asymptotic behavior is well known, see for instance the classical 
monograph by Szego, [31, Chapter VIII] or the more recent one of Ismail 
[18, Chapter 4]. Other extensions of this framework studied in the literature 
include a weight with an algebraic singularity and a discontinuity inside the 
interval of orthogonality: 

w{x) = (1 — x)“(l -|- x)^h{x)\xQ — x\'^r.c{x), 

where a,^,^ > —1, xq G (—1,1) and Sc(x) is a step function, see [10, 9]. 
In this case a separate asymptotic analysis is needed in a neighborhood 
of X = xq. It is also possible to consider Jacobi polynomials with non¬ 
standard parameters a and [3, see for instance [20, 21], but in this article we 
will restrict ourselves to the classical case, with weight function (1.1). 

2 Asymptotic expansions for Jacobi—type polyno¬ 
mials 

In this section we present large n asymptotic expansions for the orthogonal 
polynomials and related quantities. We distinguish between several regions 
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Figure 1: Regions of the complex plane in which the polynomials have 
different asymptotic expansions: the lens (I), the outer region (II) and the 
right and left disks (IH and IV). 


in the complex plane, shown in Figure 1, since the polynomial Tin exhibits 
different asymptotic behaviour in C: 

• a complex neighbourhood of the interval (—1,1) excluding the end¬ 
points, subsequently called the ‘lens’ (region I) 

• two disks around the endpoints ±1, called the right and left disk (re¬ 
gions in and EZ) 

• and the remainder of the complex plane, the ‘outer region’ (region H). 

Remark 2.1. Mathematically, the regions are of arbitrary size and, depend¬ 
ing on how they are chosen, any given point z € C can in principle belong to 
almost any region. In terms of implementation, this choice can be relevant 
as the expansion in one region may be more accurate than that in another 
region for a given point. We return to this remark in ^5.4. 

All our results are formulated in terms of a particular complex matrix¬ 
valued function R{z) G In this section we hrst elaborate briefly on 

the properties of the function R{z). Next, we introduce auxiliary functions 
that are needed in the statements of the expansions. Finally, the asymptotic 
expansions of the polynomials are stated region by region. 
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Figure 2: The system of contours T,r consists of the boundaries of the 
regions in Figure 1. The contours are oriented as shown. 

2.1 The function R{z) in the complex plane 

The function R{z) is a 2 x 2 matrix complex-valued function, that satisfies 
the following properties: 

1. R{z) is analytic (entrywise) in C\S/j, where the contour S/j is depicted 
in Figure 2. This contour consists of the boundaries of the regions in 
Figure 1. In each piece of S/j minus the self intersection points, the 
function R{z) admits boundary values R±{z), where the plus (minus) 
sign corresponds to the left (right) side with the given orientation. 

2. As n —)• oo, the function R{z) admits an asymptotic expansion of the 
form 

R{z) ^ 1 + ^ n->oo (2.1) 

k=l ^ 

which is valid uniformly for z G C \ {dUs U dUs). Here, Us and Us are 
the right and left disks respectively. 

3. The coefficients Rk{z) in the previous expansion are analytic functions 
of z in C \ {dUs U dUs). 

4. Rk{z) = 0{1/z) for z —)■ oo. 

It is crucial to note that the coefficients Rk{z) depend on 2 : and are given 
by different expressions inside and outside of the disks Us and Us- We will 
write R^^^^^{z) and R}^^^{z) to refer to the coefficients for 2 ; in the interior 
of Us and Us respectively, and R‘^'^^^'^{z) to indicate the coefficients for 2 ; 
outside these two disks. 
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Because of the above properties, this matrix R{z) is close to the identity 
as n ^ oo, uniformly in z. Thus, for the leading order behaviour of the 
expansions, one may simply substitute R{z) = I, the 2x2 identity ma¬ 
trix. Higher-order expansions are obtained automatically by determining 
the coefficients Rk{z) in formula ( 2 . 1 ) for k > \ and then plugging in an 
asymptotic expansion for R{z) that will be derived in §3 and §4. The first 
four terms are given explicitly in §A. 


2.2 Auxiliary functions 

In order to formulate the asymptotic expansions in the different regions 
of the complex plane explained before, we need some auxiliary functions. 
In this section we state their definitions. Additional comments about the 
computation of these functions are given in §5. 

The global behavior of iTniz), away from the interval [—1,1] is governed 
by the Szego function D{z) corresponding to the weight function w{z), which 
is an analytic function for z G C \ [—1,1]. In our case, because of the form 
of w{x), see ( 1 . 1 ), we can write 


D{z) 


(z-l)“/2(z + l)/^/2 /(^2 _ 1 ) 1/2 . log/l(C) dC \ 

(^( z )(«+«/2 ((2 _ 1 ) 1/2 j 


( 2 . 2 ) 


where 7 is a closed contour in the complex plane that encircles the interval 
[—1,1] once in the positive direction but not the point z, see [22, §1.1]. In 
this neighbourhood, h needs to have a positive real part and we take the 
branch of the logarithm that is real on [— 1 , 1 ]. 

We also use the function 


ip{z) = z + {z‘^ - 


(2.3) 


which is a conformal map from C\ [—1,1] onto the complex plane outside the 
unit circle. Note that we choose the branch cut of the square root on [—1, Ij. 
An alternative expression for (p{z) can be given in terms of the arccosine 
function, using the standard definition with a cut on (— 00 , — 1] U [1, 00 ), see 
[7, §4.23.24 & §4.23.25]: 

Vj(z) =exp(/0 (z)arccos(z)), 0(z) = |^’ ^ n’ (2-^) 

[ — 1 , arg(z — 1 ) < 0 . 

The function 9{z) corresponds to sgn(Imz) in C \ M, and on the real axis it 
is equal to —sgn(z — 1). This function primarily serves a practical purpose. 
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namely to ease the implementation of the branch cuts of several functions 
in this paper. 

We can rewrite D[z) as follows, using standard branch cuts: 

D{z) = w{zY^‘^ eyiY){—iO{z)'il:{z)). 


Here, the phase function 'il){z) captures the oscillatory behavior of the poly¬ 
nomials TTniz) on the interval [—1,1]. In the complex plane, it is given by 


ip{z) 


a +13 an 

-arccos z- 

2 2 




(1 - r log/r(C) dC 

Ani - 1 ) 1/2 - z' 


(2.5) 


where 7 does encircle the point z now as well, see [23, (3.9)]. The expansions 
in the next section are formulated in terms of 'ip{z), rather than in terms 
of the Szegd function D{z) itself. Note that ip{z) depends on the analytic 
function h{z) in the generalized Jacobi weight function through the contour 
integral in (2.5). 

Observe that this same contour integral. 


m{z) 


1 / log/r(C) dC 

27rz (C^ - 1)1/2 C - 2 ;’ 


is an analytic function of the variable z in C \ 7, in particular at 2: = ±1. 
Therefore, we can expand it in power series 


m{z) ^'^Cn{zm{z) ^'^dn{z + lY, ( 2 . 6 ) 

n=0 n=0 


and apply Cauchy’s integral formula to obtain 

= — I log/?-(C) dC 

27ri j; (C2-1)1/2 (C-l)"+'’ 

W log/r(C) dC ^ ’ 

27ri j; (C 2 - 1 ) 1/2 (C+ !)-+!’ 

for n > 0. These coefficients Cn and dn were introduced in [22, Lemma 6.4 
& 6 . 6 ] and are used to construct R{z), see §3.2. The convergence properties 
of ( 2 . 6 ) depend naturally on the behavior of the function h in the complex 
plane. We note the following symmetry relation: if h{—Q = h{C), then 
dn = i-ir+^Cn. 

Finally, we will need the limit of the Szegd function D{z): 

= , 14 m D(z) ^ 2 - ^ exp (T ^ ((/f ■ (2.8) 
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All contour integrals in this section can be computed either symbolically 
using residue calculus or numerically using trapezoidal rules in the complex 
plane, as explained in §5.2. 


2.3 Asymptotics of monic orthogonal polynomials Tiniz) 

2.3.1 Monic polynomials in the lens I 

Putting together the consecutive transformations in [22] for z G I, we obtain 


TTniz) = 


2 l/ 2 -r 


R 


outer 


(^) 


Doo cos(A„p(z)) 
-i/D 

OO cos (A„, 2 ( 2 ;)) 


(2.9) 


y/w{z){l — \ 0 y 

with branch cuts implemented as in §5.1 and the following phase functions: 


^n,i{z) = (n + 1 / 2 ) arccosz + ipiz) — vr/l, ( 2 - 10 ) 

^n, 2 (z) = (n — 1/2) arccos z + ipiz) — 7 r/ 4 , 

with ip{z) given by (2.5), and D^o as in (2.8). 

In particular, this expansion shows that the orthogonal polynomial has 
cosine-like behaviour in the interior of the interval (— 1 , 1 ), with a frequency 
that depends on n. One can show that the expression is actually valid in all 
of region I by analytic continuation of all underlying functions. 


2.3.2 Monic polynomials in the onter region 11 

For z G n, the asymptotic expansion is 

_ _ 2-V^-^ /1\^ ( F>ooexp(i0(z)A„,i(z)) A 

V^(l - W \-^/Doo exp (ze(z)KMz))) ■ 

( 2 . 11 ) 

Note that this formulation differs from [22, §1.2 & (9-2)]; it is numerically 
more stable, since it avoids raising ip{z) to some power and allows re-use of 
the same contour integrals. 

However, when 7 would contain points where the analytic continuation 
of log/i(z) is not guaranteed, one could use the formulas in [22, §1.2 & (9.2)]. 
It is also possible to define Xn,i{z) and \n,2(z) with the only difference that 
the contour integral in ^(z) is taken only around [- 1 , 1 ] (not z), and use 

(()%°-"(.) / „^exp(»WA..,,W) \ 

(z - l)"/ 2 (z -F l)/3/2(l - z 2 )V 4 ^i 6 »(z)An, 2 (z)^ j 



The polynomials behave like complex exponentials in region II. Note that 
away from the interval (2.11) is exponentially close to (2.9) as n —)• oo. One 
may think of the polynomials as being asymptotically a sum of two complex 
exponentials. In region I the exponentials are of comparable size and they 
combine into a cosine. In region H, one of the exponentials dominates the 
other. Hence, the asymptotic expression simplifies. 


2.3.3 Monic polynomials in the right disk IH 

A formula for x G (1 — 5,1] is given in [22, §10] and [23, (2.27)]. One can 
extend this result to z G HI: 


T^niz) = 


y/ 7 :n arccos z /1 




( 2 . 12 ) 


with 


_ /Hoo (cos (Ci( 2 :)) Ja(narccos 2 ;) + sin(Ci( 2 ;)) Jo(narccosz))\ 

^ ' (cos (C 2 (^)) arccos z) + sin ((( 2 (-s)) arccos 2 ;))y ’ 


and 


A / \ \ “TT 1 

0,2(2:) = ip { z ) + ^ ± 2 


where the + sign corresponds to Ci(2:) and the — sign to C 2 ( 2 ). 

The polynomials behave like a Bessel function near the right endpoint 
X = 1. This is typical asymptotic behaviour near a so-called ‘hard edge’, in 
the language of random matrix theory. The order of the Bessel function a 
corresponds to the order of the algebraic singularity of the weight function 
through the factor (1 — x)". It is unaffected by the presence of the analytic 
factor h{x). 


2.3.4 Monic polynomials in the left disk 12 

For 2 : G EZ, which includes the left part of the interval, we obtain similarly 

v/7rnarccos(-2) 

) — 


yO, 


R^^\z)B{z) (2.13) 


with 


B{z) = 


Doo (sm{fj,i{z))j 0 (narccos(— 2 )) -|- cos{fii{z)) (narccos(— 2 )) 

(sin(// 2 ( 2 )) J /3 (n arccos(— 2 )) + cos{fi 2 {z))J'i^ {n aiccos{—z)) 
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and 


/Ui, 2 (^) = ^ ^ ^arccosz, 

The polynomials behave like a Bessel function of order /3 near the left 
endpoint x = —1. Note that, compared to (2.12), the roles of a and j3 are 
interchanged and other symmetries can be identihed. We found that it is 
simpler to construct explicit formulas for the left disk, rather than to infer 
them from the formulas of the right disk by invoking symmetry. 


2.4 Asymptotics of leading order coefficients 

The asymptotic expansion of the leading coefficients 7 n of the orthonormal 
polynomials Pn{x) = 'Jn'^nix) is [22, §9.2] 


~ 

in 


')2n 




1 + 2iD, 


OO 

ooE 

k=l 


UlT + U'g 


(n + 1)^ 


2 , 1 / 


(2.14) 


The quantities defined and extensively described in §3. They 

are the constant 2 x 2 matrices that multiply {z^f in the expansion 

for R{z), of which we use the lower left elements here. Explicit expressions 
for these matrices up to A: = 4 are given in Appendix A. 


Remark 2.2. We do not state asymptotic expansions for the orthonormal 
polynomials. These can he obtained simply by multiplying the expansion for 
the monic polynomial with that of the leading order coefficient 7 ^. Common 
factors can be cancelled to avoid roundoff or overflow and this is included in 
the implementation. 


2.5 Asymptotics of recnrrence coefficients 

In the three term recurrence relation (1.2), the recurrence coefficients have 
the following large n asymptotic expansion [22, §9.3] 


Oin ~ 


OO 


E 


Cf + vS 

(n + 1 )^ 


+ 


1,1 


7:f‘ + 



and 


fln 


2iDl 


+E 


Cf + C'w* 


k=l 


W' 


fright ^ ^left 



1 , 2 / 
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The quantities ^ in these expressions are the same as those appear¬ 
ing in (2.14) above. Following [22, Theorem 1.10 &: §9.3], we note that the 
order 1/n terms in the previous expressions cancel out. This can be easily 
checked with the formulas for given in the Appendix A, and gives 

the estimations 

an = 0{\lr?), /3n = ^ + ^(l/n^), noo. 

2.6 Remarks on asymptotic expansions 

The asymptotic expansions are stated in [22] for points x G M on the inter¬ 
val. Proofs for the validity of their extension to points z G C in a region 
containing (part of) the interval, as they are stated in this paper, are omit¬ 
ted for the sake of brevity. One has to carefully consider the branch cuts 
involved, which are discussed in §5. 

For general ol and /3, these asymptotic expansions correspond to a relative 
error of size where T is the number of terms. If = 1/4 = /3^, all 

higher-order terms are zero {R{z) = I) and (2.9), (2.12) and (2.13) coincide 
[22, Rem. 1.14]. In that case, the leading order term of (2.9) gives already 
exponential accuracy [22, Rem. 1.5 &: 1.11] and the function h(x) is taken 
into account only in the definition of 'ip{x). If in addition h{x) = 1, then we 
obtain the explicit form of the Chebyshev polynomials. 

Although technical, these expansions can readily be differentiated and 
this is included in the implementation. In [32], derivatives were used as part 
of the computation of the points and weights of (generalized) Gauss-Hermite 
quadrature on the real line. There, in the generalized case, the polynomi¬ 
als were evaluated by a numerical solution of the corresponding Riemann- 
Hilbert problem. As we mentioned in the introduction, the expansions in 
this paper may be used to compute Gaussian quadrature rules on [—1,1]. 
In the implementation we have included a test script for this computation, 
based on a Newton method similar to that of [32]. This paper affirmatively 
answers the question raised in the conclusions of [32], whether high-order 
asymptotic expansions can be effectively derived from a Riemann-Hilbert 
formulation. An extension to Laguerre weights and generalized Laguerre 
weights is under current investigation. 

Finally, we note that the asymptotic expansions (2.9) - (2.13) also lend 
themselves to a cosine transform in order to improve accuracy near the 
endpoints (see, e.g., [3]). Accuracy may be lost in expressions involving 
(1 — x^) when x is close to ±1, due to cancellation. One may substitute 
arccosx = 6 , and then we have for example that (1 — x^)“^/^ = (sin0)“^/^. 
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which is numerically stable for 9 close to 0. In our implementation, we have 
also included series expansions around the endpoints. There, the particular 
singularity (1 — is cancelled analytically with other terms, as well as 

the singularity that arises from w{z)~^/‘^. 


3 Computation of higher-order terms 

3.1 Local jumps for the matrix R{z) 

It follows from (2.1) that the matrix-valued function R{z) is close to the 
identity matrix as n —)• oo. In fact, the leading order terms of the expansions 
in §2 are obtained by simply substituting I in the previous expressions. 

The jumps of the matrix R{z) across the contour S/j, shown in Figure 
2, tend to the identity matrix as n —)• oo, but we have two different types of 
behaviour. The first type of jump is exponentially small in n, 

R+{z) = R-{z){l + 0{e-^^^)), c>0, 

and holds on the lips of the lens-shaped region, which is the boundary 
between the regions I and H. On the other hand, on the boundary of the 
disks around the endpoints we have 

R{z) = Tl + O (3.1) 


uniformly for z G dU^ U dlls- 

The main idea to obtain higher-order terms in the asymptotic expansion 
for TTn{z) is to compute the higher-order terms Rk{z) in (2.1). To this end, 
we write the jump matrix for R{z) as a perturbation of the identity matrix, 
I -|- A(z), i.e. we write (3.1) as 

R{z) = (3.2) 


We then consider a full asymptotic expansion in powers of 1/n for A(z): 


k=l 


^k{z) 

h ’ 
r) ^ 


n OO, 


uniformly for z G The terms Ak{z) are identically 0 in T,ji\{dUsUdUs), 
because the jump of the first type is exponentially close to the identity there. 
On the boundary of the disks, the terms Afc(z) can be written explicitly as 
^nght/ieft^^^^ as we detail next. 
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3.2 The definition of Ai 


An explicit expression for Ak{z) is known [22]: 

(3.3) 

{2log cp{z)) 

X (+ I - i) - s) \ 

Here, we have used the notation (a, m) to denote, for m > 0, 


(a, m) = 


n:r=i(4«^-(2n-i)^) 

22 ™m! 


along with (a, 0 ) = 1 . Also, D^o is given by ( 2 . 8 ), and 

^ ^ “ ^/2(z2 _ 1)1/4 ^(z) V2 ) ’ 

^l( 7(z)+7(2:)-^ -i(7(2;)-7(2;)-i)\ 

2 V*(7(^;) - 7(^;)"^) 7(^;) + 7 ( 2 :)“^ y’ 


with 7 (z) = and (p{z) given by (2.3). The function F”s'^*(z) is 

F"gi^*(z) =exp [ie{z) [i^{z) + ^)) . (3.( 

This function is analytic in \ (1 — (5,1], and it has an expansion there; 


n=0 


Cn{z - l)’^ , z G Us, (3.7) 


with coefficients Cn defined in §2.2 by (2.7). In the implementation, see §4.2, 
we use this formula combined with ^{z) in terms of the arccosine, see (2.4). 
We have also used the standard notation for the third Pauli matrix (T 3 , 


1 0 
0 -1 


and fiz)"^ = 


f{z) 0 \ 

0 ’ 


for any function f{z) 7 ^ 0 . 

There are analogous formulas for A|f^*(z) for z near — 1: 
{21og|-45(z)])'= 


13 



X 


with 




i{k-k) \ 


^left(^)-<73^(^)-l^-a3^ 

(3.9) 


^left(^) 


= exp 





(3.10) 


which is an analytic function in Us \ [—1, —1 + <5), with 

{z^-l)y^J2dn{z + ir), zeUs 

n=0 / 

(3.11) 

and coefficients dn given by (2.7). 

It is important to note that by {z^ — 1)^/^ we mean the analytic branch of 
the square root that behaves like z as z —)■ oo in any direction. We comment 
further on the correct implementation of this expression in §5.1. 




exp 


Remark 3.1. The special case with = (3“^ = 1/4, mentioned before in 
^2.6, follows from the form of the matrices in (3.3) and the fact that (a, m) 
and {f3,m) vanish for any m > 1. This is easily seen from the definition 
(3.4) of (a,m). It follows also that (a,m) vanishes whenever a is a half¬ 
integer, once m surpasses a eertain maximal value. This implies that in 
these eases most eoefficients vanish identieally. This simplifies computa¬ 
tions somewhat, but not greatly, sinee R{z) does not necessarily have a finite 
number of terms nor poles. Sueh eases include Gegenbauer polynomials or 
the closely related spherieal polynomials, for example, of the kind employed 
in the spectral method presented in [26]. 


3.3 Recursive computation of Rk{z) 

We recall expression (2.1) for convenience; 

R(z) I + f" n^oo, z gC\(UsGUs). 

k=l 

The function R{z) is analytic in the regions I, IT, HI and EZ, but has jumps 
across the contour S/j. Recall that we write to refer to the 

coefficients in the interior of the right/left disks, and R^'^^^^{z) for the coef¬ 
ficients outside of the disks. 
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By expanding the jump relation (3.2) and collecting the terms with equal 
order in n, we obtain a link between the terms Rk{z) in the expansion (2.1) 
and the A^. For /c > 1, we have 


k 

^ ^ ^ fright/leftleft 

/=1 


z G dUsUdUs, 


(3.12) 


with = I, cf. [22, (8.12)]. 

This is an additive Riemann-Hilbert problem for the Rk{z)- We are 
looking for a solution to this problem recursively for each value of k, and 
at each step of the recursion, the additive jump involves the solutions of 
previous problems, i.e. Rj with j < k. 

All quantities involved are meromorphic functions in z. It should be 
noted that the functions A^.^®*^*( 2 :) and Rj{z) may have poles at z = 1, but 


^right(^) since the latter is analytic in the right disk. Thus, one 

can solve the additive Riemann-Hilbert problem as follows: 


• Expand the sum in (3.12) in a Laurent series around z = ±1. 

• Define R^“*®’’(z) as the sum of all the terms containing strictly negative 
powers of z =F 1- Since Rk{z) = 0(l/z) as z —)• oo, positive powers of 
z =F 1 do not contribute to ii^'^*®''(z). 

• Define i?^^®*^*(z) as the remainder after subtracting those poles. 

This construction ensures that is analytic outside the disk, is 

analytic inside and (3.12) holds, as required. 

A useful piece of information is conveyed by [22, Lemma 8.2]; for any 
/c > 1, the functions A^^®*^*'^^'^^*(z) have a pole at z = ±1 of order at most 
[{k + l)/2j = Ik/l], Thus, we may write 


OO 

Anght/left(z) ^ ^ =F ir, (3.13) 

m=—\k/2] 


with coefficients that can be computed explicitly by expanding 

(3.3) around z = ±1. It follows that the Laurent expansion of the sum in 
(3.12) has a principal part of the same order, which we may write as 


k 


j:^right / left ^^^^right / left 

1=1 


(z) = 


rfc/21 

E 

m=l 


j T-right / left 
^ k,m 

(z^ 1)™ 


+ 0 ( 1 ), 


z ^ ±1. 
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This expansion defines the ^ coefficients that appear in the asymp¬ 
totic expansions of the orthogonal polynomials. 

Remark 3.2. We have deseribed the coefficients at a sightly greater 

level of generality than in [22]. In the notation of [22], the first coefficients 
are 


^(1) ^ ^(1) ^ 


^(2) = _g(2) ^ 

= Ulf. 


The construction outlined above yields 


rfe/21 

= E 

m=l 


r/-right rrleft \ 

^k,m \ 

(z-l)^ (z + 1)™ j ’ 


z €C\{UsUUs). (3.14) 


At the same time, since is analytic in Us (respectively Us), it 

has a local power series expansion 


OO 

jjright / left / , 

u) ~ 2 ^ Qk,n 


(^Tl)", 


n=0 


(3.15) 


with some coefficients that can be determined as well. The three 

sets of coefficients are necessarily related. It follows from the additive jump 
relation (3.12), by expanding around 2 : = ±1 and comparing equal powers, 
that they satisfy the identities: 


fe-l ri/2l-m 

T^right / left _ -r ^right / left ^right / left-r ^right / left 

^k,m ^k^—m ' / ^ ^ '^k—j,l j^—m—l 

. /rfc /21 

=ffi I E (-^ - n + l)n(±2)--t/J;f 

fc-i ri/2i-i-n 

x^right/left ^right / left x .right / left 

“ ^k,n ~ ^k-j,l ^j,n-l ’ 

J=1 1=0 


(3.16) 

(3.17) 


where +2 corresponds to and —2 to Observe that the roles of 

the right and left superscripts are sometimes interchanged. Also, in the last 
expressions we have used the notation 


{n)m = n{n + 1) ■ ■ ■ {n + m - 1) 
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to denote the Pochhammer symbol. 

A possible approach to compute higher-order terms is to implement 
a symbolic computation of power series expansions around z = ±1, and 
to combine (3.12), (3.13) and (3.15) in order to obtain the coefficients 
■ However, this procedure turns out to be extremely inefficient 
symbolically for high-order terms, because many lengthy expressions are 
constructed and manipulated. Using the relationships (3.16) and (3.17) im¬ 
proves this situation, but in the following section we explore an alternative 
way to compute the matrices directly and more efficiently. 

4 Simplifications and explicit formulas 

The crucial formula in §3 is the jump relation (3.12), from which can be 
determined recursively. In this section, we rewrite the jump relation as (4.1) 
below, in such a way that the computation of higher-order terms is signifi¬ 
cantly accelerated. We also establish explicit formulas for the expansions of 
all quantities involved, such that higher-order terms can be computed fully 
numerically, without having to resort to a symbolic computation package. 

4.1 Simplifications 

In the computations outlined in the previous section, some combinations of 
Afc( 2 ;)’s simplify or cancel. This can be used to speed up the computation 
of considerably. We start by writing the jump relation (3.12) 

using the coefficients Rk-m{z) instead of {z). 

Proposition 4.1. The jump relation (3.12) can be written as follows: 

^right/left^^^ = Rr^^iz) - (4.1) 

m=l 

with = I and with 


m—1 

^nght/left(^) ^ ^ (4.2) 

1=1 

Proof. To prove the result, we proceed by induction for the case of the right 
disk. The base case A: = 1 is trivial and we assume that (4.1) holds until 
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k — 1. Since A: — m</c — 1, we can substitute the right hand side of (4.1) 
for in (3.12). This yields 

k / k—m \ 

m=l \ n=l / 

k k k—m 

= E -EE ^jr-“-,.(z)»?“(^)AS'*(z). 

m=l m=l n=l 

The second sum can be rewritten using a change of variables i = m + n: 

k k—m 

m=l n=l 

= E «rr(^) E »?“(^)A2“(^) = E ^"'(^) (-sTl^) + aT‘(^)) . 

£=2 n=l e=l 

where it was possible to add the case £ = 1 because s”®'^*( 2 ;) — A’(^®^*(z) = 0. 
This proves the result, and the left case is analogous. □ 


At first sight, (4.1) is merely rewriting (3.12), but this formulation has 
two essential advantages: 

• The jump term in (4.1) is written in terms of rather than 

and the former has a simple and non-recursive expression (3.14). 

• The definition of the coefficients ^ can be greatly simplified to 
a non-recursive expression too, involving just the A^’s. 

More precisely, we have the following result: 

Proposition 4.2. The terms defined by (4.2) satisfy 

fright/left(_^) ^ 


for odd m and 




left 

^m 



4a^ -f 2m — 1 (a, m — 1) ^ 
’ \n{ip{z))^ 2™+im ' 

4/3^ -|- 2m — 1 (/3, m — 1) ^ 
ln(-v9(z))™ 2"^+im 


for even m, with {a,m) defined by (3.4). 


This can be proven again by mathematical induction, see the proof in 
Appendix B. 
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4.2 Precomputing a series expansion for 

The recursive procedure in §3.3 relied on series expansions. Symbolic ma¬ 
nipulation of series, which might not be available or slow, can be avoided by 
deriving explicit formulas for the expansion. In this section, we construct 
the expansion for the functions in order to find the 

In view of Proposition 4.2, this is equivalent to deriving the expansion of 
the Afc’s. In turn, this amounts to deriving expansions for all quantities ap¬ 
pearing in their (lengthy) definition in §3.2 and combining these expansions 
through convolutions to obtain the final result. This is conceptually straight¬ 
forward, but laborious in practice. In this section, we supply rather many 
technical details, as great care has to be taken with signs and branch cuts, 
which can only be achieved with thorough understanding of the methodology 
of [22]. 

We want to compute the coefficients II4,m in 

OO 

E z^±l. (4.3) 

m=—\k/2\ 


We proceed by detailing the expansion of the quantities appearing in 
definition (3.3) and afterwards, one by one. We observe that one can write 
(3.3) as 


fright/left 


(g,fc - 1) 

(21og[±^(z)])^ 




(4.4) 


where q = a for the right disk and q = j3 for the left disk. Here, the ± signs 
always correspond to the right/left endpoint. The function Gk{z) in (4.4) 
can be given in terms of M{z), see (3.5), and see (3.6) and 

(3.10). Omitting superscripts for brevity, we have 


Gkiz) = M{z)F{zr «) F{z)-^^M-\z) (4.5) 

with 

Working out the multiplication of the matrices for odd and even k, we obtain 




1 

(z2 - 1)1/2 


—az 

ia 



+ ,b( “4!/a+«) 

'^-zcos(y„+^_i) 


-icos(2/«+/3+i)\ 

-cos(y„+^) J 
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and 


^even/ \ = (^ 0\ b ( -sin(2/a+^) ism{ya+0+i)\ 

^ Vo oy (2:2 - 1)1/2 Visin(y„+^_i) sin(y„+^) )' 

Based on (3.7) and (3.11), the functions = y^{z) above are given by 


2/7 ~ log(±</3(2)) - i{z^ 


00 


1)'/^ E 

n=0 


Cn{z - 1)” 

dn{z + l)^, 


(4.6) 


with 7 = Q; + /3or7 = a + /3±l. 

In order to compute the coefficients Wk^m in (4.3), we will expand all the 
previous functions in power series around 2 = ±1. We will use the notation 
u = 2 2F 1, with minus (plus) sign for the right (left) disk. 

We start with the power log(±(/7(2))“^ in (4.4); from (2.3), we get 


log (p{z) = ^0(2) arccos(2), 


(4.7) 


and log(—(/9(2)) = \ogip{z) — 6{z)TTi. Expanding the arccosine as 2 —)■ ±1, 
we obtain 

00 .lx 

log(±^7))~{±27'-'^E/„i.”, /•. - ■ ns) 

using the standard Pochhammer symbol {^)n and the variable v explained 
before. We note that the factor 9{z) in (4.7) is cancelled by the branches of 
the logarithm and the square root. Continuing, we have the recursive result 


(log(±V?(2))) ^ ~ (±2u) 

n=0 


with 5:1^0 = I//0 = 1) and, for k > 1, 


9i,n — 


-1 

~h 


n—1 

j=0 


(log(±(^(2)))-'= ~ (±2 u)-'=/2 

n=0 


9k,n — ^ '^ 9k—l,l9l,n—l- (4.9) 

1=0 


In order to expand cos(|/xy) and sin(y.y), we note first that because of 
(4.6) and (4.8), we have 


2/7 ~ -/(± 2 u)^/^ ^ Pl,n,'yV'^, 

n=0 


Pi,12,7 



Cn-j2 1 

d„_7-2)-i 
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Note that with the standard branch cuts for the powers, is real on 
the interval [—1,1]. Then, for k > 1, 

oo n 

V'y ^ ^ ^ ^ Pk^n^'-y'^ •> Pk^n,'y ^ ^ P/c—^,7; 


n=0 


and 


cos 


n=l 


sm 


n=0 r 

00 

Vi ~ -i{±2vf/‘^ ^ 


^ m 


1=0 


V^, 


eveuyn = 


E(±2) 

j=0 


j /^2j+l,n-j,7 


(2j + 1)! 




d(±2n)i/2^ 

?i=0 n=0 

One more expansion is needed, as we have to divide by — 1)^/^. Since 

oo 

n 

n=U ' 

we obtain 


(^2 _ i)-i/2 ^ (±2r;)-V2 (±2) 

n=0 ^ 2^ 


cos(y. 


7^ 


(z2 - 1)1/2 

sm(j/7) 
(z2 - 1)1/2 


(±2r;)-i/2^ 


n=0 


n . 1 


i=i 


1 + E / (±2^^; 


.7 zjodd 


V 


— I 


>E 


n=0 

Also, we observe that 
z 


” /-i\ 


i=o 


E( /j(±2)-^i^; 


— 7 Tj-even 


72 _ 1)1/2 


(±2r;)“^/2^('2n + 1) 


n- 2 

n 




n=0 

to complete the computation of G^'^'^(2;) and 

Finally, bearing in mind (4.4) and Proposition 4.2, we write the coeffi¬ 
cients Wk^m as follows. 

Proposition 4.3. The coefficients *11 expansion (4.3) for the 

functions are given explicitly by 


W, 


m+{k+l)/2 

right/left _ W; ^ Ij ,^odd 


Km (±2)3^2 


E! 9k,jG 


k,m+{kd-l)/2—j > 


j=0 
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fright/left _ 1) I {4:Q -\-2k ^)9k,m^k/2 


k^m 


(± 2 ) 3^2 


2k 


m-\-k/2 

-^+ ^ 9k,j 
j=0 


^even 

^k,m+k/2—j 


with q = a for the right disk and q = /3 for the left disk. Here, g^j are 
defined by (4.9) as the coefficients in the expansion o/log(±(/?(z))“^ around 
z = HI. The coefficient matrices are the expansion coefficients of 

Gk{z) defined by (4.5) around ±1 for odd and even k, respectively. 

Remark 4.4. One can compute the values directly in a simi¬ 

lar way: compare (3.13) and (4.3), using the correspondence between these 
expressions given in Proposition 4-2. 

Analogous formulas to (3.16)-(3.17) can be derived, relating the 

and values. However, the coefficients can 

also be used to compute directly, based on (4.1). This requires 

fewer values than (3.16)-(3.17) uses values. That leads 

us to the final formula: 


k-l \{k-j)/2\ 
i=l i=max(m—l'j/2] ,1) 


k-l \j/2\-m 

+E E 

j=l n=0 


T(fc-1)/21 

E 

2 = 1 


(1 i ?T-)n Tj-left / right 

O i 


(± 2 )* 


k-j,i 


Ti^right / left 
j,—n—m 

{±2Yn\ 


5 Numerical issues and implementation 

5.1 Square roots and other algebraic singularities 

Several multivalued functions appear in the asymptotic expansions of §2, 
whose implementation in the complex plane deserves some attention. Recall 
first the mathematical expression for the ip function, first introduced in (2.3), 
which is 

ip{z) = z + {z‘^ - 1)^/^. 

This function is understood to be analytic in C \ [—1,1] and to behave 
like z as z —7- oo. This means that (z^ — ig the analytic continuation 
of the square root Vx^ — 1, positive for x > 1, to the complex plane minus 
the interval [—1,1]. Observe that the square root is negative when z < — 1, 
on the negative real axis. 
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This poses a problem in implementation, since the standard branch cut 
of the square root function in (z^ — 1 )^/^ results in an extra cut on the 
imaginary axis, because the argument of the square root is real and negative 
there. This extra cut is avoided when implementing the expression 

phi(z) = z+sqrt(z- 1 )*sqrt(z+ 1 ) 
using the standard branch cuts. 

Similar considerations apply to other multivalued functions such as 
(1-z2)1/2^ (1-z 2)V4 and _ l)V4^ 


which are understood as analytic continuation of the corresponding functions 
on the real axis. 

Finally, the arccosine function appears repeatedly in §2, including in 
expression (2.4) for (p(z), dehnition (2.5) of tp(z), dehnition (2.10) of A±(z), 
and in the expansions for the polynomials. The standard arccosine function 
has a branch cut on (—oo,—1] U [l,oo). On the cut, the boundary values 
are the following, see [7, 4.23.24 & 4.23.25]; 


=Filog((z^ - + z), 2 :G[ 1 ,oo), 

TT =F i log((z^ — 1 )^/^ — y, z £ (—oo, — 1 ]. 


5.2 Computation of contour integrals 

Several expressions in §2 involve contour integrals around the interval [—1,1], 
see (2.8), (2.5) or (2.7). The general form of these integrals is 

where 7 encircles the interval [— 1 , 1 ] once in the positive direction and is 
contained in the region where h{z) is analytic and has a positive real part, 
see [ 22 , § 1 . 1 ]. 

If log/i(C) appearing in F’(C) is an entire function, or meromorphic with 
known poles, these integrals can be computed explicitly using residue cal¬ 
culus. For instance, if log/i(C) is entire, we only need to pick up the residue 
at infinity: 

±.Jnc)<ic = -F.u 
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where -F_i comes from the Taylor-Laurent series expansion 


oo 

~ Fmt^, t —7- 0. 

m=—OO 

An interesting example for which the coefficients Cn and dn can be com¬ 
puted in this way occurs for the weight function h{z) = exp(—with 
m > 1, see §6.1. 

If this approach is not possible due to lack of analyticity of F{(), these ex¬ 
pressions can be evaluated with the trapezoidal rule along a suitably chosen 
contour. This technique is exponentially accurate in the number of function 
evaluations, since the integrands are analytic and periodic functions, see 
[33]. We propose to integrate along Bernstein ellipses; 

T'p = { I 0 G [0, 27r] I , P > 1 

These are parameterized by a value p > 1, with p = 1 corresponding to the 
interval [—1,1] itself and p > 1 to an ellipse with foci at ±1. The size of the 
parameter p is limited by the analyticity of the integrand in the complex 
neighbourhood of [—1,1]. It may be possible to determine an optimal value 
of p: we refer to [4] for an extensive analysis of the optimal radius in circular 
Cauchy integrals and to [35] for a related study of the optimal value of p of 
Bernstein ellipses in the computation of Chebyshev coefficients. For (2.5), 
the countour also has to encircle the point x at which we wish to evaluate 
V’(x). 

An explicit expression for the trapezoidal rule using M points is 
/ F{QdC, = [ F' (sPe*® + ^0 

J E p J 0 

k=0 

with equispaced points located at 9^ = 2'Kk/M. The minimal number of 
points M to use is of course dependent on the integrand. Due to the expo¬ 
nential convergence of trapezoidal rules for periodic integrands, the number 
M can in general be taken to be fairly small, except in the vicinity of poles 
of the integrand. The successive doubling algorithm in [4] that gives an 
optimal M (which should increase with n in Cn and dn) is included in the 
implementation. 
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We note that care has to be taken in general to remain on the same 
branch of the analytic continuation of log h{z). In other words, if we have 
Im(log/i(C)) ^ (—vr, vr], evaluating h{z) first and then taking the principal 
branch of the logarithm would not yield the correct answer. In that case, 
one could take p closer to 1 and a higher M, or better, fill in the analytical 
continuation of \ogh{z) into the trapezoidal rules. 

5.3 On the analyticity and positivity of h{z) 

Several contour integrals in §2 are given in terms of the logarithm of h, so 
it is instructive to understand its possible behaviour in the complex plane. 
Recall that the principal branch of the logarithmic function has a branch 
cut along the negative real axis. 

The function h satisfies several conditions stated in [22]: 

1 . /i is a real-valued and positive function on [—1,1], 

2 . h is analytic in a complex neighbourhood of [—1,1], 

3. furthermore, the real part of h is strictly positive in a complex neigh¬ 
bourhood U of [—1,1]. The contours in §2 are restricted to lie in U. 

The first condition guarantees existence of the orthogonal polynomials for 
all n. The second condition is required for the complex deformations in the 
Riemann-Hilbert problem to be valid. 

We elaborate on the third condition. First, if h vanishes at a point 
on [—1,1], then the asymptotic behaviour of the orthogonal polynomials 
becomes substantially different. For examples of such behaviour, see e.g. 
[10, 9]. Second, if the real part of h has positive and negative values in 
a region, then there may be a branch cut of the principal branch of the 
logarithm in that region. In particular, branch points arise at roots of h{z) 
in the complex plane. Though branch cuts may be moved, and the contours 
appearing in this paper may be deformed in order to avoid branch points 
and other singularities of log h, the simplest implementation uses Bernstein 
ellipses confined to the region where h{z) has positive real part. 

It may appear to be problematic that h{z) appears in the asymptotic 
expansions of the polynomials in the complex plane through w{z)~^^'^ when 
h{z) has singularities there. Clearly, the polynomials do not have such 
singularities. However, one may verify that singularities of h{z) cancel and 
the asymptotic expansions are, in fact, analytic functions away from the 
interval: see §2.3.2 for example. 
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5.4 Sizes of the region 

To conclude, we return here to Remark 2.1 about the sizes of the different 
regions of the complex plane. Since the sizes of the disks around the end¬ 
points and the size of the lens can be chosen arbitrarily, different expansions 
can be valid at any given point in the complex plane. We observed ex¬ 
perimentally from our heuristics test in the implementation that whenever 
different expansions are valid at a point, the corresponding relative errors 
in the approximation of the polynomial differ typically at most by a factor 
of about 2, for large n. 

There are a few exceptions. First, the expansion in the outer region 
is less accurate as we approach the interval [—1,1] since it is ‘missing’ one 
of the exponentials that combines into the cosine-like expression in (2.9). 
This ‘missing exponential’ is exponentially small in the outer region, hence 
it can be discarded there, but not inside the lens. The point x = 0.2 -|- 0.5z 
is not on the interval, but close to it, and indeed in Figure 4 we only see 
the expected (order of) accuracy for the expansion in the outside region 
starting from n = 32 for the highest number of terms in Figure 4. The 
exponentially small difference between (2.9) and (2.11) is only negligible 
from there onwards. 

Another exception to the factor 2 difference appears when evaluating in¬ 
side a disk of radius about 0.2 around the endpoints. There, the expansions 
in the respective disks can be orders of magnitude more accurate than the 
ones in the other regions. Indeed, the latter expansions blow up at the end¬ 
points, whereas we note that in (2.12) and (2.13) the singularities ai z = ±1 
are only apparent, something that is reflected in the series implementation. 
We also remark that when evaluating very close to an endpoint, say at a dis¬ 
tance where Cm is the machine epsilon, one needs the series expansion 
(3.15) of (which also avoids the explicit subtraction of poles 

that happens in (4.1)), as well as a series expansion of the other factors 
in (2.12) and (2.13). Without the use of series expansions, it is certainly 
helpful to employ the cosine transform, as commented on before in §2.6. 

6 Examples and numerical results 

An important source of examples is given by the canonical modifications 
or perturbations of the Jacobi weight function, via polynomial or rational 
factors (Christoffel and Geronimus -with mass equal to 0- perturbations). 
See [18, §2.7] for further references. In this section, we illustrate the accuracy 
of the asymptotic expansions with three different examples, that were chosen 


26 



from literature. 


6.1 An exponential weight function 

Consider the weight function 

w{x) = h{x) = exp (—, 


with a = /3 = 0. This weight function appears in methods for avoiding the 
Gibbs phenomenon of Fourier series [12], 

The residue calculus from §5.2 yields explicit formulas for the coefficients 
Cn and dni because \ogh{z) = —cz^'^ is an entire function. We have that 
Cn = dn = 0 for n > 2m — 1. An explicit formula for the other coefficients is 


Cn 


—C 



2m — 1 — 2j 
2m — n — 1 — 



By symmetry, we have that dn = (—l)’^“'"^Cn. 
Also, 


V’(x) = - ^Q!(arccos x — tt) + (3 arccos + 




X" 


2m—1 

E 

n=0 


Cn(x - 1)*' 


and 


= 2 -“/ 2 -/ 3/2 


exp 


—cfm — 1/2 


m 


We illustrate the accuracy of the asymptotic expansion in the left disk. 
Figure 3 shows the relative error at the point x = —0.97 as a function of 
n, for c = 7 and m = 2. The ‘exact’ polynomials we compared with were 
obtained using Matlab routines from the OPQ-library that accompanies the 
book [11]. It is clear from the hgure that the expansions improve with 
increasing n, at a rate that depends on the number of terms. High accuracy 
is achieved already at moderate values of n, for example 10“^ relative error 
is seen at n = 32 using six or seven terms. However, for small n, expansions 
with fewer terms are more accurate than expansions with more terms, as is to 
be expected from the asymptotic nature of the expansions. The asymptotic 
expansions of the coefficients 7n, On and f3n, in the other regions and for 
other values of x exhibit similar behaviour. 
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Figure 3: Relative error of the asymptotic expansion in the left boundary 
region as a fnnction of n, for the weight function w{x) = exp (—The 
expansion is evaluated at x = —0.97, with a varying number of terms. For 
each number of terms i, a line is plotted with slope that interpolates 
the best relative error. 


6.2 A Jacobi-type weight function with a branch point in 
the complex plane 

Next, we consider the weight function 

w(x) = - -- - , 

a/(1 - x)(x + 3) 

which leads to a = —1/2 and /3 = 0. It appears in the approximation of 
non-periodic functions on an interval using Fourier series on a larger interval 
[17, §3]. 

In this case, log/i(z) is not entire dne to the singularity at x = —3. We 
have used the trapezoidal rules explained in §5.2 in order to compute the 
relevant contour integrals. We chose a Bernstein ellipse with p = 4, which 
crosses the real axis at x = —2; that is halfway between the singularity 
X = — 1 of the integrand in and related quantities and the singularity at 
—3. This choice reduces roundoff errors, although computing the condition 
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number like in [4] and [35] seems to give an optimal p very close to 6 as 
predicted there. It would suffice to use only M = 80 points in the trapezoidal 
rule, which is in between the last two iterations M = 64 and 128 computed 
here by successive doubling up to C 2 - 




Figure 4: Relative error of the asymptotic expansion in the lens (left) and 
outside region (right) as a function of the degree n, for the weight function 
w{x) = l/v^ — x)(a: + 3). Both expansions are evaluated at the same 
point X = 0.2 + 0.5z and for a varying number of terms. 

Figure 4 shows that we still obtain the expected order of convergence of 
the relative error. Some saturation appears for the highest number of terms 
around 10“^^, due to doing computations close to machine precision and 
accumulating errors in the recurrence relation for the ‘exact’ polynomials. 
Results are shown for the asymptotic expansion in the lens as well as in the 
outer region, but with both expansions evaluated at the same point x = 
0.2 + 0.5i. Such comparisons may lead to a decision as to which expansion 
to use in which part of the complex plane, see §5.4. 

6.3 Toda measures 

Our results include the Toda modification explained in [18, §2.8] and given 
by h{x) = with t G M. The resulting orthogonal polynomials appear 
in the literature as time-dependent Jacobi polynomials, and they have been 
studied in connection with integrable systems and Painleve transcendents, 
see for instance [1]. 

For this weight function, we have that cq = —t = do (which only enter in 
the third terms of the expansions), c„ = 0 = for n > 1, Uoo = 
and 'i/j{x) = ^ ( a(arccos x — vr) + /3 arccos x j — |\/l — x^- The leading order 
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term of the orthonormal polynomial in the lens is 


Pn{x) ~ 


1/2 cos ^[n + (1 + a + /3)/2] arccos(x) — 7r/4 — a7r/2 — tVl — x‘^/2'j 
-y/w(l — exp(—xt/2) 



Figure 5: Leading order terms of the orthonormal polynomials in the lens 
for the weight function w{x) = (1 — exp(—xt) at n = 101. 

Figure 5 shows us that the envelope of the polynomial indeed behaves as 
The weight function will become very small at x = 1 when t —)> + 00 , 
making the polynomial ill-defined and large there, while p„(—1) will become 
very small. The inverse is true for t —>• —00 and the cases t = —2 and t = 2 
are symmetric. In this example we have chosen a = j3 = —1/2, and in this 
case we can simply use the expansion in the lens throughout the interval. 
The relative error with respect to the true polynomials remains bounded by 
10“^ pointwise for all cases shown. 
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A Expressions for the first four higher-order terms 


We illustrate the recursive computation of Rk by giving the first few terms 
explicitly for 2 : outside the two disks. We have: 
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and 


y42(a, b,c) = + 8a + 86 + 8c — 46^ + 1, 

B2{a, b,c) = — 8a — 86 — 8c + 4 a^ + 46 ^ — 10, 

6*2(0, 6, c) = — 8a — 86 — 8c — 4 a^ — 46 ^ + 10, 

D2{a, b,c) = — 8a — 86 — 8c — 46 ^ + 1 . 


Next, we have 

.rright _ 4 q;^ - 1 3 / A3{a,P,co,-do) i{q3 + r3){a, l3,co,-do)\ 

8192 \i{q3-r3){a,ld,co,-do) D 3 {a, jd, cq, -do) J °° 


^Note that the following expressions are slightly different from those given in [22, §8.2). 
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2 /4(a, b, fo, go, fi) [48(46^ - 1 ) 50(50 - /o - 3(a + b)) + b^{48a^ + 542) 

+48/0^(46^ + 120^ - 28) + 144(a + b ){ Ab ‘^ + 80^ - 19)/o + 0^(566^ + 422) 
+06(1152(0^ + b ^) + 988a6 - 2880) + 16a® - 13936^ - 9510^ - 498 + 86®] . 
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B Proof of Proposition 4.2 

The basic idea for the proof is induction, and also the fact that the term 

{sj{z)Am-j{z) + S,rr_j(2:)A4z))"®’"*/'®^* 

always simplifies to 0 when m is odd or something proportional to I when 
m is even. 

Proof. For simplicity we present the case of sf^^^{z) and omit the super¬ 
scripts. The proof proceeds by induction in m in formula (4.2). The case 
m = 1 is clear from (4.2) because then there is an empty summation in 
(4.2). We assume that the proposition holds for /c = 1,..., m — 1. 

If m is even, we can rewrite the sum on the right hand side of (4.2) as 
follows: 

m-l m/2-l 

'y ] Sj{z)A.rn—j{z) = ^ + Sm—j{z)Aj{z')] , 

i=i 3=1 
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and if m is odd, 

m-l (m-l)/2 

^ Sj{z)A^_j{z) = [Sj{z)AYn—j{z) + Sm—jiz)Aj(z)] . 

j=i i=i 

We recall the formula (3.3), and for 1 < A: < m — 1, we denote 


and 


Then 


^k = 


Bt = 




k odd, 
k even. 


2k 


Sj(^z')Am—j{z)-\-Sm—j(^z')Aj(^z') — 


{a,j - l){a,m- j - 1) 


(21ogv9(z))™ 

where = BfA^_- + B^_-A^. 

We observe that if m is odd then for j odd we have m — j even, and for 
j even we have m — j odd. In both cases, a direct computation shows that 
= 0. Then Sm(z) = Amiz), and the proposition is true in this case. 
If m is even, then j and m — j are simultaneously odd or even. In this 
case, the matrix reduces to a multiple of the identity matrix: 


where 

Therefore, 


— \ ■ T 

_^(4a^ +4j(m-j) - l)(4a^ - {2j - l){2{m-j) - 1)) 

8j{m-j) 


Sj{z)Am-j[z) + Sm-j{z)Aj[z) — (2 log (B-l) 

Now, we need to sum this last expression over j. Using the symmetry 
j ^ m — j of the coefficients, we can write 

m/2-l 

^ml 2 {Z')Arfij 2 {z) + ^ ^ [Sj(z)(z) + Sm—j{z)Aj(^z')] 
i=i 
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(B.2) 


— 2 [Sj{z)Am-j{z) + Sm-j{z)Aj{z)] 

J = 1 

^ m—1 

= 2(2log^(z)r 

using (B.l). The sum (B.2) is of hypergeometric type, since we can write 
the coefficients (a, k) as follows: 

(a, k) = (_i)fc ( 2 +«)j 2 (2 3 ) 

in terms of the standard Pochhammer symbol. Next, we may apply a known 
algorithm due to Gosper and later extended by Zeilberger, see [14] and also 
[19]. Let 

«i = («,J - 

then the algorithm seeks Sj, such that the sum telescopes as follows: 


m—1 

aj = Sm-1 - Sq. 

i=i 

Under the hypothesis that Sj/Sj-i is a rational function in j, the ratio 
ttj/aj-i is also rational in j, and can be written as 

Q-j _ Pj 

Pj-iG 

where in this case 

Pj = (4a^ — 4j^ + 4jm — l)(4a^ + 4j^ — 4jm + 2m — 1), 

Qj = {-m + j - l){2a + 2j - 3)(2a - 2j + 3), 
rj = {2a + 2{m — j) — l)(2a — 2{m — j) + l)j. 

Then Sj is constructed as follows: 

S, = (B.4) 

Pj 

where fj is a function of j to be determined. In this case, fj is a polynomial 
of degree 2 in j that satisfies the linear recursion 

Pj ~ 'L+i/i “ G'/j-i- 
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From here we get 


4a^ — + 4jm — Aj + 2m — 1 


and using (B.4), we get Sj. Finally, a brief computation yields 

Sm-l Sq - 

4a^ + 2m — 1 , , 

=-(a, m — 1), 

m 

using (B.3) and the fact that m is even. This completes the proof for 

This reasoning can be carried out analogously for the by replacing 
a by (3, ln((/9(2:)) by \n.{—ip{z)) and taking into account the extra sign of the 
off-diagonal elements in □ 
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